##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--maf_cancer_file"), type = "character") ,
    make_option(c("--maf_im_file"), type = "character") ,
    make_option(c("--maf_cancer_msi_file"), type = "character") ,
    make_option(c("--maf_im_msi_file"), type = "character") ,
    make_option(c("--images_path"), type = "character") ,
    make_option(c("--info_file"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    maf_cancer_file <- paste(work_dir,"/maf_public/All_use.maf",sep="")
    maf_cancer_msi_file <- paste(work_dir,"/maf_public/All_use.msi.maf",sep="")
    maf_im_file <- paste(work_dir,"/maf_public/All_use.IM.maf",sep="")
    maf_im_msi_file <- paste(work_dir,"/maf_public/All_use.IM.msi.maf",sep="")
    images_path <- paste(work_dir,"/images/mutRate",sep="")
    info_file <- paste(work_dir,"/public_ref/combine/MutationInfo.combine.addMolecularSubType.tsv",sep="")


}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

maf_cancer_file <- opt$maf_cancer_file
maf_im_file <- opt$maf_im_file
info_file <- opt$info_file
images_path <- opt$images_path
maf_cancer_msi_file <- opt$maf_cancer_msi_file
maf_im_msi_file <- opt$maf_im_msi_file

dir.create(images_path , recursive = T)

###########################################################################################

info <- data.frame(fread(info_file))
dat_maf_cancer <- data.frame(fread( maf_cancer_file ))
dat_maf_im <- data.frame(fread( maf_im_file ))
dat_maf_cancer_msi <- data.frame(fread( maf_cancer_msi_file ))
dat_maf_im_msi <- data.frame(fread( maf_im_msi_file ))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
## 注释肿瘤样本的类型
maf <- dat_maf_cancer
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_use <- merge( maf_use , info[,c("Tumor" , "Class" , "Molecular.subtype")] , by = "Tumor")
maf_cancer <- maf_use

maf <- dat_maf_cancer_msi
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_use <- merge( maf_use , info[,c("Tumor" , "Class" , "Molecular.subtype")] , by = "Tumor")
maf_cancer_msi <- maf_use

###########################################################################################
## 肠化样本
maf <- dat_maf_im
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_use <- merge( maf_use , info[,c("Tumor" , "Class" , "Molecular.subtype")] , by = "Tumor")
maf_im <- maf_use
maf_im$Class <- "IM"

maf <- dat_maf_im_msi
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_use <- merge( maf_use , info[,c("Tumor" , "Class" , "Molecular.subtype")] , by = "Tumor")
maf_im_msi <- maf_use
maf_im_msi$Class <- "IM"

###########################################################################################
## 总的，不同来源的
## IGC和DGC分开的，合并的
maf_cancer2 <- rbind( maf_cancer , maf_cancer_msi)
maf_cancer2$Class <- "GC"
maf_use <- rbind(maf_cancer , maf_cancer_msi , maf_cancer2 )
maf_use <- rbind(maf_use , maf_cancer , maf_cancer_msi , maf_cancer2 , maf_im , maf_im_msi)

###########################################################################################
## 提取recurrent位点
maf_use$vid <- paste(maf_use$Hugo_Symbol , 
    maf_use$Chromosome , maf_use$Start_Position , 
    maf_use$Reference_Allele , maf_use$Tumor_Seq_Allele2 , sep = ":")

###########################################################################################

info_gc_divide <- subset( info , Class != "IGC + DGC" ) 
info_gc <- info_gc_divide
info_gc$Class <- "GC"
info_im <- subset( info , From == "NJMU" & Class != "IGC + DGC" ) 
info_im$Class <- "IM"

info_use <- rbind( info_gc_divide , info_gc , info_im )
info_all <- info_use
info_all$From <- 'All'
info_use <- rbind(info_all , info_use)

###########################################################################################
## 按照病理类型计算突变率
## 以人为单位
maf_use$Molecular.subtype <- ifelse( maf_use$Molecular.subtype %in% c("EBV" , "unknown") , "EBV/unknown" , maf_use$Molecular.subtype )
maf_use$Molecular.subtype <- ifelse( maf_use$Molecular.subtype %in% c("MSI" , "POLE") , "MSI/POLE" , maf_use$Molecular.subtype )
info_use$Molecular.subtype <- ifelse( info_use$Molecular.subtype %in% c("EBV" , "unknown") , "EBV/unknown" , info_use$Molecular.subtype )
info_use$Molecular.subtype <- ifelse( info_use$Molecular.subtype %in% c("MSI" , "POLE") , "MSI/POLE" , info_use$Molecular.subtype )

mut_rate <- maf_use %>%
group_by( Hugo_Symbol , vid , Class , Molecular.subtype) %>%
summarize( MutNum = length(unique(Tumor )) )
mut_rate$id <- paste0( mut_rate$Class , ":" , mut_rate$Molecular.subtype )

#mut_rate <- subset(mut_rate , MutNum > 1)
class_num <- info_use %>%
group_by( Class , Molecular.subtype) %>%
summarize( SampleNum = length(unique(Tumor)) )
class_num$id <- paste0( class_num$Class , ":" , class_num$Molecular.subtype )

mut_rate <- merge( mut_rate , class_num[,c("id" , "SampleNum")] , by = "id" )
## 分IGC、DGC和IM
mut_rate$MutRate <- mut_rate$MutNum/mut_rate$SampleNum

mut_rate_final <- mut_rate

###########################################################################################
## 至少2人的IM或GC
record <- subset( mut_rate_final , Class %in% c("IGC" , "DGC" , "IM"))
record <- subset( record , MutNum > 1 )

result <- subset( mut_rate_final , vid %in% record$vid  )

out_name <- paste0( images_path , "/MutRate.RecurrentPoint.molType.tsv" )
write.table( result , out_name , row.names = F , quote = F , sep = "\t" )